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I.  INTRODUCTION 

This  report  addresses  the  problem  of  generating  a  surface-fitted  grid  in 
a  model  fin-body  problem  consisting  of  a  circular  cylinder  with  four  identical 
symmetric  fins  attached.  This  grid  is  to  be  used  in  the  calculation  of 
incompressible,  laminar  flow  at  moderate-to-high  Reynolds  numbers.  The  aim  of 
the  calculation  is  to  resolve  the  details  of  the  separated  zone  at  the  leading 
edge  of  the  fin-cylinder  juncture  and  the  subsequent  vortex  that  forms 
downstream.  Thus,  the  grid  must  have  proper  clustering  so  as  to  resolve  the 
regions  of  high  flow  gradients. 

The  approach  used  here  is  to  generate  the  grid  analytically  but  to 
determine  the  metric  coefficients  numerically.  Such  an  approach  has  been 
pursued  successfully  by  Jameson  [1]  and  Caughey  and  Jameson  [2-4]  in  solving 
three-dimensional  inviscid  transonic  flows  about  wing-body  combinations.  The 
basic  idea  is  to  map  the  physical  geometry  to  a  strip  of  almost  constant  width 
using  a  sequence  of  conformal  transformations.  Then  boundary  fitted 
coordinates  are  generated  by  the  application  of  a  shearing  transformation. 

The  result  of  the  latter  transformation  is  a  nonorthogonal  coordinate  system 
in  the  physical  plane  but  one  in  which  the  nonorthogonality  can  be  controlled. 

The  present  work  is  an  extension  of  the  Jameson-Caughey  technique  for 
what  is  called  the  wind  tunnel  problem  to  the  case  of  an  initial  value  plane 
ahead  of  the  airfoil.  In  order  to  treat  viscous  flow,  clustering  trans¬ 
formations  are  used  so  that  the  computational  grid  is  uniform  in  all  three 
directions . 

One  advantage  of  the  present  technique  is  that,  owing  to  the  simple 
cylindrical  body  geometry,  a  three-dimensional  grid  is  generated  by  stacking  a 
series  of  two-dimensional  grids.  Another  advantage  of  the  analytical  approach 
over  the  numerical  solution  of  elliptic  partial  differential  equations  as  a 
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means  of  grid  generation  is  its  much  greater  speed  which  is  especially 
important  for  three-dimensional  applications. 

II.  ANALYSIS 

2.1  Geometry  of  Computational  Domain 

We  start  the  grid  generation  analysis  by  defining  the  geometry  about 
which  a  surface  fitted  grid  is  to  be  generated  and  the  extent  of  the 
computational  domain. 

1.  The  body  is  an  infinitely  long,  hollow  circular  cylinder  of 
radius  Rc  with  its  centerline  parallel  to  the  free-stream 
velocity  vector. 

2.  Four  identical  fins  of  constant  unit  chord  and  infinite  span, 
consisting  of  symmetric  airfoil  sections,  are  mounted  on  the 
cylinder  90  degrees  apart  with  their  chord  planes  passing 
through  the  cylinder  axis. 

3.  The  computational  domain  consists  of  the  region  interior  to 
an  outer  cylinder  of  radius  Rt  which  encases  the  inner 
cylinder  and  fins,  bounded  upstream  and  downstream  by 
planes  normal  to  the  cylinder  axis. 

A  schematic  of  one  fourth  of  the  geometry  and  computational  domain  is  shown  in 
Fig.  1  and  a  head-on  view  showing  the  coordinate  system  in  the  crossflow  plane 
appears  in  Fig.  2.  Since  the  fins  are  identical  and  equally  spaced,  we  have 
four  planes  of  symmetry,  namely,  at  9  3  0,  n/4,  tt/2  and  3ti/4.  Thus,  in  the 
flow  field  calculation  for  this  model  problem  only  the  segment  0  <  9  <  it/4 
needs  to  be  considered. 
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2.2  Sequence  of  Transformations 

Four  transformations  applied  In  sequence  are  required  to  map  the  fin- 
cylinder  and  surrounding  computational  domain  into  a  rectangular 
parallelpiped.  Then  a  fifth  stretching  transformation  is  applied  to  adjust 
the  grid  line  spacings  for  proper  flow  field  resolution  in  physical  space  and 
to  allow  a  uniform  step  size  in  all  three  computational  coordinates. 

We  start  by  defining  polar  coordinates  (r,9)  in  the  crossflow  plane,  as 
shown  in  Fig.  2,  according  to 

r  =  (y2  +  z2)1/2  ,  (1) 

0  =  tan-*  .  (2) 

Thus,  points  in  physical  space  are  defined  by  standard  cylindrical  coordinates 
(x, r, 0 ) . 

Following  Caughey  and  Jameson  [2] ,  the  first  transformation  normalizes 
(x,r,0)  according  to  (all  lengths  are  referred  to  the  airfoil  chord): 


x  =  x  -  ds  +  2n  2  ,  (3) 


0  =  49  ,  (5) 

where  ds  is  the  location  of  the  singular  point  of  the  unwrapping  trans¬ 
formation  and  is  just  inside  the  leading  edge  of  the  airfoil.  Note  that 
in  the  above  definitions,  0  <  r  <  1  and  0  <  9  <  Tt  in  the  computational 
domain.  The  upper  limit  on  9  is  convenient  in  the  next  transformation. 


****** 
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Because  r  =  constant  is  a  surface  fitted  coordinate  we  need  only  generate 
a  surface  fitted  grid  in  the  (x,9)  plane.  The  geometry  of  an  r  =  constant 
surface  in  the  computational  domain  is  sketched  in  Fig.  3. 

The  conformal  transformation 

x  -  16  =  £n[l  -  cosh(£  +  ip)]  ,  (6) 


applied  to  an  r  =  constant  surface  unwraps  the  geometry  in  Fig.  3  to  produce 
the  domain  shown  in  Fig.  4.  The  minus  sign  has  been  used  on  the  left  in 
Eq.  (6)  so  that  the  upper  symmetry  plane  maps  to  the  positive  £  axis. 

In  the  present  problem  initial  conditions  from  an  axis  etric  boundary 
layer-potential  flow  composite  solution  are  specified  on  tl  'lane  x  =  -  a. 
This  initial  value  line  in  an  r  =  constant  surface  (IVL)  is  Dwn  as  line 
segment  ABC  in  Fig.  3.  Under  transformation  (6),  the  IVL  it  -  .o  a  near  semi¬ 
circle  in  the  (£,n)  plane,  as  shown  in  Fig.  4.  The  airfoil  image  in  this 
plane  is  the  arc  DEF . 

We  next  apply  another  conformal  transformation  to  nearly  straighten  out 
the  IVL  in  Fig.  4.  This  transformation  is 


£  + 


in  =  £  +  in  + 


So2 

£  +  in  ’ 


(7) 


where  £Q  is  the  intersection  of  the  IVL  with  the  £  axis  (Point  A  in  Fig.  4). 
The  conformal  transformation  (7)  maps  the  upper  and  lower  boundaries  in  the 
(£,n)  plane  into  slowly  varying  functions  of  £  in  the  (£,n)  plane,  as  shown 
in  Fig.  5.  We  note  that  near  Points  A  and  C  the  IVL  is  now  cusp-like. 
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The  fourth  transformation  Is  a  shearing  transformation  which  straightens 
out  the  upper  and  lower  boundaries  in  the  (£,n)  plane.  This  transformation 
i  s 


X  =  K  . 


(8) 


Y 


(9) 


7  =  r  ,  (10) 

where 


D  =  D(C,r)  =  nlT  -  nL 


(11) 


and  Pp  aod  are  the  ordinates,  at  a  given  £,  of  the  upper  and  lower 
boundaries  in  the  (S,n)  plane. 

Finally,  to  provide  for  clustering  the  grid  lines  near  the  fin  and 
cylinder  surfaces  to  resolve  the  viscous  layers  there  and  to  space  lines 
around  the  airfoil  and  in  the  wake  as  desired,  we  introduce  one-dimensional 
stretching  functions  as  follows: 


=  Fj(x;  , 

(12) 

=  f2(y)  , 

(13) 

=  f3(z)  , 

(14) 

For  the  time  being  we  leave  Fj,  F2  and  F-j  unspecified.  Thus  (Xc,  Y  ,  7-c )  are 
the  computational  coordinates  devised  so  that  the  step  sizes  AXC,  AYC  and 
A 7,c  are  constants. 
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2.3  Conformal  Mapping  Relations 


Since  the  FORTRAN  code  Is  written  in  terms  of  real  variables,  the  real 
and  Imaginary  parts  of  the  conformal  mappings  must  be  determined.  In 
addition,  the  inverses  of  both  mappings  are  needed  because  the  grid  generation 
procedure  requires  being  able  to  proceed  from  the  (x,9)  plane  to  the  (Xc,  Y  ) 
plane  and  then  back  to  the  (x,0)  plane. 

The  real  and  imaginary  parts  of  Eq .  (6)  yield  the  two  relations: 


cosh  E,  cos  n  =  1  -  ex  cos  0  , 


sinh  £  sin  n  =  ex  sin  6  .  (If 

The  solutions  for  x  and  0  are  obtained  by  squaring  (15)  and  (16),  then  adding 
and  making  use  of  the  ordinary  and  hyperbolic  trigonometric  identities.  The 
result  for  x,  choosing  the  proper  sign,  is 

x  =  £n(cosh  £  -  cos  n)  ,  (17 

and  0  is  obtained  from  Eq.  (15),  viz., 

0  =  coS-l  fi.  -  cosh  g  cos  n1 

cosh  £  -  cos  n  v 

To  obtain  the  solutions  for  £  and  p  we  first  define, 

p  =  1  -  ex  cos  9  ,  (19 


q  =  ex  sin  9  .  (2C 

Following  the  same  procedures  as  above,  we  eliminate  n  to  obtain  a  quadratic 


ca^t^r 
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equation  for  sinh^  £  which  has  the  solution 


sinh2  f  =  j  [  (  S2  +  4q2)l/2 


-  B] 


(21) 


where 


B 


1  "2  "2 

1  -  p  -  q 


(22) 


In  the  right  half  plane  £  is  the  positive  root  of  F.q.  (21).  The  expression 
for  n  with  the  proper  behavior  (0  <  n  <  it)  is  obtained  from  Eq.  (15),  viz. 


n  =  cos 


-1 


f  _ 

p 

1 

cosh 

j 

(23) 


Next,  the  real  and  Imaginary  parts  of  Eq.  (7)  yield 

C  2 


£  =  £ 


n  =  n 


1  + 


1  - 


£2  +  n2  ’ 


Co 


2  ) 


£2  +  n2  ’ 


(24) 


(25) 


We  determine  £Q  from  Eq.  (17)  by  setting  x  =  -  a  =  -  a  +  In  2  and  n  -  0.  The 
The  result  is 

£0  =  cosh” ^ ( 1  +  2e~a)  .  (26) 


where  a  =  dg  +  dIVL* 

To  solve  for  £  and  n  in  terms  £  and  n  we  return  to  the  complex  form  which 
is  written  as. 


w  =  z 


(27) 


where 
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w  =  \  +  in  , 

(28) 

z  =  C  +  in  . 

(29) 

Solving  Eq.  (27)  for  z  yields, 

T  i  9  2 

r  =  4  w2  -  (0  » 

(30) 

where 

1 

♦  =  z  -  j  w  . 

(31) 

Let  us  now  define 

4>  =  u  +  lv  . 

(32) 

Then,  combining  Eqs.  (28),  (29)  and  (31)  gives: 

e  =  u  +  j  l  , 

(33) 

n  =  v  +  y  n  . 

(34) 

Now  Eq.  (30)  leads  to  the  following  relations: 

u2  -  v2  =  p  , 

(35) 

uv  *  q  , 

(36) 

where 

1  ,-2  -2.  2 

P  =  U  -  n  )  -  , 

(37) 

q  =  4  ^  • 

(38) 

Equations  (35)  and  (36)  can  be  solved  for  u  and  v  with  the  result: 
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u2  =  j  (p  +  p)  ,  (39) 

v2  =  ~  (p  -  p)  ,  (40) 

where 

P  =  (p2  +  hq1)112  .  (41) 

Then  the  final  result  for  £  and  n,  combining  Eqs.  (33),  (34),  (39)  and  (40), 
is 

C  -  Y  l  +  (U  +  p)]1/2  ,  (42) 

n  =  y  n  +  ty  -  p)ll//2  •  (43) 

2 .4  Calculation  of  Shearing  Boundaries 

The  shearing  boundaries,  which  are  straightened  out  hy  the  shearing 
transformation  Eq.  (9),  are  defined  as  9^(1)  and  n^(C)*  Thus  is  the  image 
of  the  upper  airfoil  surface  and  the  line  9=0  downstream  of  the  trailing 
edge  while  9l  is  the  image  of  the  upper  half  of  the  initial  value  line 
(x  =  -  a)  and  the  line  9  =  it  for  x  >  -  a. 

We  start  by  determining  the  image  of  the  upper  half  of  the  airfoil  in  the 
(S,p)  plane.  The  airfoil  will  be  given  as  a  set  of  points  (xp,  yp)p  where  for 
convenience  we  take  the  origin  at  the  leading  edge.  Then  the  scaled  airfoil 
coordinates  in  the  (x,9)  plane,  for  a  given  r,  are: 

xp  =  Xp  +  Jin  2  -  dg  ,  (44) 
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8p  =  4  sin-^ 


Next,  the  image  in  the  (£,n)  plane  is  computed  from 


y? 


(45) 


where 


£p  =  sinh"1  [I  (a  -  8) ] 1 /2 


hp 


cos 


-1 


cosh  £p 


(46) 

(47) 


a  =  (82  +  4q2)1/2  ,  (48) 

and  p,  q  and  0  are  given  by  Eqs.  (19),  (20)  and  (22).  Then  the  image  in  the 
(£,n)  plane  is 


£f  =  ?{•(!  +  * 

(49) 

nu  =  hp(l  -  u)  , 

(50) 

and 


u 


Co" 


?T 


(51) 


The  upper  boundary  beyond  the  airfoil  trailing  edge  is  the  image  of  0  =  0 
which  maps  to  n  *  n.  To  calculate  riy  in  this  region  we  first  compute  a 
uniform  point  distribution  of  £  on  the  interval  (CTE'Cmax)»  Then 


! 
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corresponding  values  of  £  are  computed  by  iteration  from 


r  (n+1 )  =  _ § _ 

1  +  P*  ’ 

where  superscript  n  denotes  the  iteration  number,  and 


(52) 


We  note  that  Eq.  (52)  converges  quite  rapidly.  With  a  value  of  £  known, 
is  computed  from 


(53) 

r'U 


nu  =  ir*(l-p).  (54) 

In  the  calculation  of  the  lower  n  boundary  the  shearing  transformation 
requires  that  the  same  £  distribution  be  used  as  was  determined  for  ny.  The 
lower  boundary  is  computed  in  two  segments,  the  first  on  the  Interval  (0,£o), 
where  £0  is  the  image  of  £0,  and  the  second  on  the  remaining  interval 
^o’^max^ 

On  the  interval  (0,£o)  we  calculate  £  and  n  by  iteration  from  the 
rapidly  convergent  formula: 


f(n+l  )  _  _ § _ 

*  1  +  u 


(55) 


where  in  this  case 


y 


So" 


2  x  ( n) 


(£"  +  n4) 


(56) 


n  =  cos~^(cosh  £^n)  -  2e_a) 


(57) 


To  start  the  Iteration  we  set  u  =  1  in  Eq.  (55)  which  from  Eq .  (56)  is  seen 
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to  be  exact  at  £  =  £Q.  With  £  and  n  known,  nL  is  calculated  from 

nL  =  n  •  (1  “  P)  •  (58) 

On  the  interval  (£0,  £max)  we  know  from  Eq .  (58)  that  the  image  of  6  =  tt 
is 


nL  =  0  .  (59) 

Thus  knowing  the  distribution  of  nu  and  on  (0»!max)  we  can  obtain  the 
distribution  of  the  shearing  distance  D  from  Eq.  (11). 

2.5  Stretching  Functions 

The  approach  taken  here,  as  already  mentioned,  is  to  use  one-dimensional 
stretching  functions,  as  indicated  by  Eqs.  (12),  (13)  and  (14).  In  the  present 
application  the  location  and  length  scales  of  regions  of  rapid  variation  of  the 
solution  are  known,  beforehand.  In  a  Z  =  constant  plane  of  the  computational 
domain,  as  shown  in  Fig.  6,  clustering  of  Y  =  constant  lines  is  needed  near 
Y  =  1  and  0  to  resolve  the  boundary  layer  developing  on  the  airfoil  and  the 
region  around  the  corner  singularity,  x=-a,9=it,in  the  physical  plane. 
Thus,  for  the  variable  Y  a  two-sided  stretching  function  is  required.  Because 
of  the  primary  viscous  layer  on  the  cylinder  clustering  is  needed  near  Z  =  0 
which  requires  a  one-sided  stretching  function  for  Z.  The  stretching  function 
for  X  depends  on  criteria  related  to  the  flow  field  and  the  mapping  geometry 
which  will  be  discussed  later. 

Vinokur  (5)  has  determined  approximate  criteria  for  the  development  of 
one-  and  two-sided  stretching  functions  of  one  variable  which  give  a  uniform 
truncation  error  independent  of  the  governing  differential  equation  or 
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difference  algorithm.  He  Investigates  several  analytic  functions  but  finds 
that  only  tan  z,  where  z  is  real  or  pure  imaginary,  satisfies  all  of  his 
criteria. 

We  start  with  the  stretching  function  for  Y  and  note  that  both  Y  and  Yc 
are  normalized  variables  as  required  in  Vinokur's  functions.  In  the  present 
case,  z  is  taken  to  be  pure  imaginary  which  leads  to 

tanh(YcA<t>) 

^  A  sinh  A$  +  (1  -  A  cosh  A4>)  tanh(YcA(J>)  ,  (60) 

where 

A  =  (SQ/Si)172  ,  (61) 


J/2 


B  =  (SqSj ) ' 

and  Sq  and  are  dimensionless  slopes  defined  as 


s0 


dY< 

dY 


(0)  , 


(62) 


Si 


> 


which  control  the  clustering  at  Y  =  0  and  Y  =  1,  and  A<j>  is  the  solution  of  the 
following  transcendental  equation: 


R 


sinh  Aft 
A$ 


(63) 


To  avoid  solving  Eq.  (63)  by  iteration,  Vinokur  determines  the  following 
extremely  accurate  approximate  solutions  for  small  and  large  B: 


For  B  <  2.7829681 
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A4>  =  (6B)1/2(1  -  0.15B  +  0.057321A29B2 

-  0.024907295B3  +  0.0077424461B4 

-  0.0010794223B5)  , 


where 


B  =  B  -  1  . 


For  B  >  2.7829681 


(64) 


(65) 


A<J>  =  V  +  (1  +  l/VHn(2V)  -  0.02041793 


+  0.24902722W  +  1.9496443W2  -  2.6294547W3 
+  8.5679591 1W^  ,  (66) 


where 


V  =  In  B  , 


(67) 


and 


W  =  1/B  -  0.028527431  . 


(68) 


An  example  of  this  two-sided  stretching  function  for  Sq  =  100  and  =  10  is 
shown  in  Fig.  7.  For  this  case,  A<J>  computed  from  Eq.  (66)  is  5.926. 

The  one-sided  counterpart  of  Eq.  (6)  is  antisymmetric  about  the  mid-point 
and,  in  terms  of  Z  and  Zc,  is  given  by 


Z 


1  + 


tanh  [|-  AiJ>(Zc  -  1) ) 
tanh 

2 


0  <  Z  <  1 


(69) 


s--  -Vr,  - 
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where  A$  Is  the  solution  of 


SO  = 


sinh  A » 
A  <j> 


and 


S0  = 


(0)  . 


(70) 


Two  examples  of  this  one-sided  stretching  function,  Sq  =  10  and  100,  are  shown 
in  Fig.  8. 

The  stretching  function  in  x  is  required  to  have  the  following 
properties: 

(1)  It  must  have  the  ability  to  cluster  points  near  the  nose  of  the 
airfoil  to  resolve  rapid  flow  field  variations  in  that  region. 

(2)  Control  points,  where  grid  lines  are  required,  are  the  corner, 

X  =  XQ,  and  the  airfoil  trailing  edge,  X  =  Xt£. 

(3)  Downstream  of  the  airfoil  trailing  edge  where  flow  gradients 
are  decreasing  the  step  size  should  gradually  increase. 

(4)  The  stretching  function  should  have  continuous  first 
derivatives . 

(5)  For  proper  flow  field  resolution,  the  number  of  steps  on  the 
intervals  (0,XQ)  and  (XQ,XTE)  are  to  be  parameters. 

The  above  requirements  dictate  the  stretching  function  be  made  up  of  three 
piecewise  continuous  segments  on  (0,XQ),  on  (x0»XTg)  and  on  (XTE ,Xmax) . 

We  start  by  defining  variables  normalized  on  the  corner  location, 


An  appropriate  stretching  function  on  the  first  segment  is  given  by  Eq.  (61)  of 
Vinokur,  viz. 
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X  =  Xc  (1  +  j  (S0  -  1X1  -  Xc)(2  -  Xc))  ,  0  <  Xc  <  1  , 


where  Sq  is  the  slope  at  the  origin  and  is  used  to  control  clustering  of  points 
in  that  region.  The  uniform  step  size  on  Segment  1  is  given  by 

-  sr  ■  (73) 

where  is  the  number  of  intervals  on  Segment  1.  We  note  that  AXC,  as  given 
by  Eq.  (73),  is  also  the  step  size  on  Segments  2  and  3. 

On  Segment  2,  the  scaled  trailing  edge  coordinate  is  given  by 

(XC)TE  =  1  +  N2  AXc  ,  (74) 

where  N2  is  the  number  of  intervals  on  Segment  2.  We  note  that  (XC)^E  *  X,^ . 

The  constraints  to  be  satisfied  by  the  stretching  function  of  Segment  2  are: 

...  1 

X  =  1,  X*  =  Xj  on  xc  =  1 


X  =  XXE  on  Xc  =  (Xc) 


where 


y  »  =  dX 

A.  ~ 

dxc  Xc=l 


which  from  Eq.  (72)  is 


Xj  =  jO  -  s0)  .  ( 

With  three  constraints  a  parabola  is  appropriate.  The  resulting  stretching 
function  is 


"vs  - 
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X  =  1  +  [X*  +  A(XC  -  1 ) ) ( X c  -  1)  , 


a: 


where 


A  = 


XTE  -  1  -  X{[(XC)TE  -  1] 


t(xc)TE  -  1 


(7! 


On  Segment  3  a  geometric  progression  is  used  to  increase  the  step  size 
in  X.  Requiring  continuity  of  X  at  the  junction  with  Segment  2,  we  have 


xk  ~  XTE  +  AX1 


1  -  C 


k-1 


1  -  C 


,  k  >  2  , 


(7' 


where  C  is  the  constant  step  size  ratio  defined  by. 


AX. 


C  =  >  1  . 

"k-1 


Continuity  of  the  first  derivative  at  the  junction  is  ensured  by  choosing  AXj 
equal  to  the  last  AX  on  Segment  2.  No  attempt  is  made  to  match  Xmax  exactly. 

The  stretching  function  for  X  is  seen  to  have  four  parameters,  Sq,  Nj , 
N2  and  C,  which  provide  considerable  flexibility  in  the  point  distribution  of 
X.  A  typical  example  is  shown  in  Fig.  9. 
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3.  RESULTS  AND  DISCUSSION 


3.1  Generation  of  the  Grid 


The  step-by-step  procedure  to  generate  a  grid  in  the  physical  plane  for 
given  airfoil  shape  and  initial  value  plane  location  is  as  follows: 

(1)  The  uniform  computational  grid  (Xc  ,YC  ,ZC  )  is  first 

i  j  k 

established  and  then  are  calculated  via  the 

stretching  functions  described  in  Section  2.5. 

(2)  With  (X^  Yj.Z^)  known,  r^  is  determined  from 


rk  =  zk 


Then  for  r  fixed,  the  points  in  the  (X,Y)  plane  are  transformed 
to  the  (C,n)  plane  by 


where 


^ijk  =  Xi  » 

nijk  =  YjDik  +  ^L^ik  * 


Dik  ~  <nU>ik  "  (nL>ik 


By  Eq .  (45),  0p  depends  on  r  and  hence  r  and  therefore  Py 
and  Py  must  be  computed  anew  for  each  value  of  r.  The 
procedure  used  here  Is  to  calculate  more  points  than 
needed  on  the  shearing  boundaries  for  a  given  r  and  then 
to  use  Lagrange  cubic  interpolation  to  determine  Py  and 
Py  for  a  given  £. 


30  March  1983 
GHH : lhm 


(3)  With  (C,h)  known,  the  transformation  to  the  (?,n)  plane  is 

,  _  ,  -  1 /2 
Sijk  "  i  ?ijk  +  Iy  ^  +  P^^  ’ 


1  -  i  -  1/2 

hijk  "  i  nijk  +  ^2  P')  va 


where 


p  =  \  (l2  ~  n2)  -  Co2  . 


q  =  r  5  h 


(4)  Next,  the  points  in  the  (£,n)  plane  are  transformed  to 


the  (x,0)  plane  by 


xijk  *  «.n(cosh  Sijk  -  cos  nijk)  , 

(  I 

_  _i  1  “  cosh  C|^k  cos  nijki 

9ijk  "  cos  cosh  Cijk  -  cos  n1jk  |  * 


(3)  The  final  step  is  to  compute  the  cylindrical  coordinates  of 
each  grid  point  from: 


xijk  xijk  +  ds  -  in  2  , 


9ijk  =  4  9ijk  > 
rk  -  Rc  +  ^t  -  ^c^rk  * 
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3 . 2  Features  of  the  Grid 

The  shearing  transformation  applied  at  the  fourth  stage  necessarily 
produces  a  nonorthogonal  grid  in  the  (x,9)  plane.  The  nonorthogonality  is 
smallest  on  the  lower  shearing  boundary,  under  most  conditions,  and  largest  at 
the  airfoil  surface  on  the  upper  shearing  boundary,  as  can  be  seen  from  Fig.  5. 
On  the  upper  (airfoil)  boundary  the  nonorthogonality  near  the  leading  edge 
(5  =  0)  can  be  controlled  by  proper  location  of  the  singularity  of  the 
unwrapping  transformation,  F.q.  (6).  Away  from  the  leading  edge  the  only  control 
over  nonorthogonality  is  to  keep  the  airfoil  reasonably  thin,  say  eight  percent 
or  less,  which  will  maintain  ri(j  as  close  to  the  image  of  n  =  tt  as  possible. 

The  parameter  which  controls  grid  orthogonality  near  the  airfoil  leading 
edge  is  dg  in  Eq.  (3).  The  most  nearly  orthogonal  system  in  this  region  is 
produced  when  the  leading  edge  maps  into  an  n  =  constant  line.  In  the  (x,B) 
plane  such  a  line  is  closely  approximated  by  a  parabola  centered  about  9=0  and 
is  effectively  characterized  by  its  radius  of  curvature  at  the  origin,  given 
by 


Po 


d^x 


d0 


9=0 


(94) 


We  determine  pu  by  setting  q  =  n^E  =  constant  in  Eqs.  (13)  and  (16), 
differentiating  the  result  twice  with  respect  to  9  to  find  d  x/d9  ,  plus  noting 
that  dx/d9  =  0  at  9=0  and  by  virtue  of  Eqs.  (3)  and  (5)  that 
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i 

! 


The  result  is 


.  ,  -  1_  Sl"2  ^ 
p°  16  cos 


(95) 


From  Eq.  (3)  evaluated  at  the  airfoil  leading  edge  (x  =  0)  we  have 


XLE  =  £n  2  ~  ds 

and  from  Eq.  (17)  with  E  =  0  and  n  =  n^E  we  find  that 


-d. 


cos  nLE  =  1  ~  2e  s 


from  which  it  follows  that 


sin  hle  =  2[e  ds(l  -  e  ds)]1^ 


Hence,  Eq .  (95)  for  p0  becomes 


po  =  g 


d  g; 


1  -  e  s 


2e  ds  -  lj 


(96) 


(97) 


(98) 


(99) 


which  can  be  solved  for  d<-  to  yield, 


ds  =  in 


r  ^ 

1  +  16  pJ 


1  +  8  Poi  ' 


(100) 


Next,  we  fit  the  airfoil  leading  edge  by  an  osculating  parabola,  viz. 


x  =  KG'1 


(101) 


where  K  =  xj_/9j.  and  (x^,0^)  are  appropriate  airfoil  coordinates  near  the 
leading  edge.  The  radius  of  curvature  of  the  airfoil  at  the  leading  edge  is, 
from  Eq.  (101 ), 


(102) 
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The  optimum  value  of  ds  (which  produces  the  most  nearly  orthogonal  grid  near 
£  =  0)  is  obtained  by  equating  pLE  and  pu.  Thus,  dg  can  then  be  determined 
from  Eq,  (100).  Figure  10  shows  the  variation  of  ny  with  £  for  a  six  percent 
thick  Joukowsky  airfoil  for  three  values  of  ds ,  one  of  which  was  determined  by 
Eqs.  (100)  and  (102).  In  these  three  cases,  we  have  dg  <<  d^  which  has  the 

effect  of  limiting  the  influence  of  dg  on  riy  to  the  region  0  <  £  <  £u  where  here 

£0  =  0.87.  As  r  increases  from  Rc  to  Rt  the  leading  edge  radius  of  curvature 
of  the  airfoil  decreases  because  9p  decreases — see  Eq .  (45).  Thus  dg  must  be 
decreased  accordingly. 

On  the  lower  shearing  boundary  the  nonorthogonality  arises  from  the 
mapping  of  the  initial  value  line  (IVL)  by  Eq.  (7).  In  the  (£,q)  plane  the  IVL 

is  very  nearly  half  of  an  ellipse  with  the  ratio  of  the  semi-major  to 

semi-minor  axes  lengths,  defined  as  X  =  q0/£0  (nQ  is  the  value  of  n  on  the  IVL 
at  (  =  0)  given  by 


\  -  cos~l(l  -  2e  a) 
cosh~*(l  +  2e-a) 


(103) 


Figure  11,  in  which  X  is  plotted  versus  ”a",  shows  that  as  "a"  becomes  large 
X  approaches  unity  and  therefore  the  IVL  approaches  a  semi-circle  in  the  (£,q) 
plane.  Thus  for  to  have  the  smallest  maximum  (at  £  *  0)  and  hence  for 
£  =  constant  lines  at  q  =  to  be  as  nearly  orthogonal  as  possible,  "a"  should 
be  large,  say  3  or  4,  a  circumstance  desirable  on  physical  grounds  anyway. 
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At  the  Image  of  the  airfoil  trailing  edge  in  the  (£,n)  plane  (Points  D 
and  E  in  Fig.  5)  when  the  trailing  edge  angle  is  finite  the  derivative  of  ny 
with  respect  to  £  wiii  be  discontinuous.  At  the  ends  of  the  IVL  (Points  A  and 
C  in  Fig.  5)  the  behavior  of  is  cusp-like  which  means  that  the  second 
derivative  of  with  respect  to  £  is  discontinuous.  These  discontinuities 
produce  similar  type  discontinuities  in  Y  =  constant  lines  via  the  shearing 
transformation.  This  behavior  is  one  of  the  disadvantages  of  algebraic 
mappings  involving  shearing  transf ormat ion  which  is  absent  in  grids  generated 
by  solving  elliptic  partial  differential  equations.  The  discontinuous  behavior 
of  derivatives  of  Y  =  constant  lines  in  the  physical  plane  should  therefore  be 
accounted  for  in  the  calculation  of  affected  metric  coefficients  and  in  the 
numerical  method  of  solution  of  the  viscous  flow  equations. 

3.3  Numerical  Examples 

For  simplicity  a  symmetric  Joukowsky  airfoil  was  used  in  the  numerical 
examples  of  the  grid.  The  ordinates  of  this  airfoil  (for  unit  chord)  are 
given  by, 

yp  =  2/7  a  -  xf)(4  xFd  -  xF)l  '  ,  (104) 

where  xp  is  measured  from  the  airfoil  leading  edge  and  t  is  the  maximum 
thickness  to  chord  ratio.  Two  example  grids  in  the  (x,9)  plane  are  presented 
with  parameters  listed  in  Table  1  below.  The  parameter  J  is  the  number  of 
points  In  the  Y  direction. 


.  /  ■ 


-28- 


30  March  1983 
GHH : lhm 


Parameter 

Case  1 

Case  2 

Ni 

15 

n2 

— 

15 

J 

31 

31 

dIVL 

3.0 

3.0 

ds 

0.05 

0.05 

d0B 

3.0 

3.0 

T 

0.12 

0.06 

Rc 

1.0 

1.0 

r/Rc 

1.0 

1.0 

s0 

10 

10 

Si 

10 

10 

s0 

— 

0.2 

c 

— 

1.2 

Table  1.  Grid  parameters  for  Numerical  Examples 

Case  1  is  shown  in  Fig.  12  and  Case  2  in  Fig.  13.  Case  1  has  no  stretching 
function  in  X  and  no  X  =  constant  line  through  the  corner.  The  non¬ 
orthogonality  of  the  grid  in  Case  1  (12%  thick)  is  seen  to  be  more  pronounced 
at  the  airfoil  surface  than  in  Case  2  (6Z  thick)  which  bears  out  the  remark 
made  earlier.  Notice  that  both  examples  are  for  the  grid  in  the  (x,0)  plane 
on  the  cylinder  surface  (r  =  Rc)  which  corresponds  to  the  intersection  of 
the  fin  with  the  cylinder.  Hence  in  these  examples,  by  Eq.  (45),  the  airfoil 
thickness  in  terms  of  0  is  a  maximum  and  thus  the  nonorthogonality  is  most 
pronounced. 

The  computer  code  listing  is  given  in  the  appendix. 


t 

1 

i 

i 
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Figure  1.  Schematic  of  Geometry  and  Computational  Domain. 
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Figure  2.  Coordinate  System  in  Crossflow  Plane. 
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Figure  9.  Three  Segment  Stretching  Function  for  X. 
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10  C  PROGRAM  NAMES  CGHID3 

20  C  THIS  PROGRAM  COMPUTES  A  SURFACE  FITTED  C-GRID  FOR  A  FIN 

30  C  CYLINDER  8Q0Y. 

40  C  THE  FIN  IN  THIS  VERSION  IS  A  SYMMETRIC  JOUKO*SKY  AIRFOIL. 

50  C  THIS  IS  THE  3-D  VERSION. 

60  c****************************************************************** 

70  IMPLICIT  REALMS  (A-H,0-Z) 

80  COMMON  /BLKOi/  IMAX, JMAX, ITE, ITEM, ILAST,ISEGl,IStG2 

90  COMMON  /BLK02/  XlfcM , XIO , X I  BO 

100  COMMON  /8LK03/  Cl , C2 , Ci , C4 , C5 , PI , PISJ 

110  COMMON  /oLKOm/  XF(101)  ,  YFUOl) 

120  COMMON  /6LK05/  XIB(151),bbARC151),ETARL(l51) 

130  COMMON  /BLKOe/  SYO , SY 1 , SZO , SXO , SSR 

140  COMMON  /BLK07/  ZC C 1 51) , bIGZ C 1 5 1 3 

150  C 

160  1  FORMATC  514  3 

170  2  FORMATC5F10.43 

180  10  F0RMAT(1H1,4X, 'INPUT  PARAMETERS  f OR  C-GR1D'3 

190  11  FORMA  T(lriO,9X,  'ISEG1  =  '  ,  I6/10X ,  '  ISEG2  *  '  ,  16/ 1  OX ,  *  JMAX 

200  1 16/ 1  OX , ' KMAX  =  ' I6/10X, 'ITE  =',163 

210  12  FOKMATCIOX, 'DIVL  = ' , Fl 0 . 4/ 10X , ' DOB  =',F10,4/ 

220  110X, 'TAU  =» ,F10.4/10X, ' RC  = ' , Fl 0 . 4/ 1  OX . ' RT  = • , FI  0 . 4/ 1  OX  , 

230  2'SYO  =' ,F10.4/10X,  'SY1  = » , Fl 0 . 4/ 1  OX , ' SZO  = ' , F 1 0 . 4/ 1  OX , ' SXO 

240  3F1U.4/10X, ' SSR  =‘,F10.4) 

250  13  FORMAT! 1H03 

260  14  F0KMATC5X, 'STACKED  C-GRID  FOR  FIN-CYLINDER  GEOMETRY' 3 


270 

13  F0RMATC10X 

DS  =',U14.4) 

280 

C 

290 

C 

INPUT  REQUIREMENTS 

300 

C 

310 

C 

ISEG1 

s 

NO.  INTERVALS  ON  FIRST  X-SEGMENT. 

320 

C 

ISEG2 

«• 

NO.  INTERVALS  ON  SECOND  X-SEGMENT. 

330 

C 

IMAX 

= 

NG.  PGINTS  IN  X-UIRECTION. 

340 

C 

JMAX 

= 

Nu.  POINTS  In  Y-DIRECTION. 

350 

C 

KMAX 

= 

NO.  POINTS  IN  Z-OIKECTION. 

360 

C 

ITE 

= 

NO.  POINTS  ON  AIRFOIL  INITIALLY, 

370 

C 

DIVL 

3 

DISTANCE  FROM  AIRFOIL  L.E,  TO  INITIAL  VALUE  LINE 

380 

C 

DS 

3 

DISTANCE  FROM  AIRFOIL  L.E.  TO  SINGULARITY  OF 

390 

C 

COORDINATE  SYSTEM, 

400 

C 

DOh 

3 

DISTANCE  FROM  AIRFOIL  f.E.  TO  OUTFLOW  BOUNDARY. 

410 

C 

TAU 

S 

AIRFOIL  MAX.  THICKNESS  TO  CHORD  RATIO. 

420 

C 

RC 

3 

INNER  CYLINDER  RADIUS,  IN  TERMS  OF  AIRFOIL  CHORD 

430 

C 

RT 

3 

OUTER  CYLINDER  RADIUS,  IN  TERMS  OF  AIRFOIL  ChOkD 

440 

C 

SYO 

3 

Y-STRF TCHING  PARAMETER  AT  AIRFOIL  SURFACE. 

450 

C 

SY1 

3 

Y-STRETChING  PARAMETER  AT  INITIAL  SURFACE. 

460 

C 

SZO 

3 

Z-STKETCn ING  PARAMETER  AT  INNER  CYLINDER. 

470 

c 

SXO 

3 

INITIAL  X-STRETCHI NG  PARAMETER,  SEGMENT  1. 

480 

c 

SSfi 

3 

X-GEOMETRIC  PROGRESSION  RATIO,  SEGMENT  3. 

490 

c 

500 

READ (5,13 

ISEG1 , ISEG2 , JMAX, KMAX, ITE 

510 

READ( 5,23 

DIVL, DOB 

520 

READ ( 5,23 

TAU,RC,PT 

530 

READ (5, 2 3 

SYO,SY1,SZO,SXO,SSR 

540 

ITfc.MsITE-1 

550 

WRITEtb, 10) 

560 

WRITE ( b , 1 1 ) 

I5EG1,ISEG2, JMAX, KMAX, ITE 

570 

RRITEC  6 , 1 2 ) 

DIVL, DOB, TAU, RC,RT, SYO, SY1, SZO, 5X0, SSR 

580 

WRITE (b, 13) 

590 

WRITE ( b , 1 4 ) 

600 

c 

610 


CJ3?.ODO*TAU/DSORT(27.0DO) 


PI*3,14159265358979D0 

PISQ=PI*PI 

XE=1.0UO+DOB 
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620 

630 

640 

650 

660 

670 

680 

690 

700 

710 

720 

730 

740 

750 

760 

770 

790 

790 

800 

810 

820 

030 

840 

050 

860 

670 

880 

690 

900 

910 

920 

930 

940 

950 

960 

970 

980 

990 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 

1100 

1110 

1120 

1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 


C 

C  CALCULATE  AIRFOIL  COORDINATES. 

C 

CALL  FOIL 
C 

C  CALCULATE  ZC  AND  BIGZ. 

C 

CALL  STRFZ(ZC,6IGZ,KMAX,SZ0) 

DELRsRT-RC 

C 

C  6EGIN  CALCULATION  OF  STACKED  GRID. 

C 

DO  50  K=1,KMAX 
RAD=RC+OELR*BIGZ(K) 

C 

C  CALCULATE  DS  -  DISTANCE  FROM  AIRFOIL  LEADING  EDGE  TO 

C  SINGULARITY  Of  UNWRAPPING  TRANSFORMATION. 

C 

THr sDASINCYFC4)/RAD) 

FHOs0.5D0*THF*THF/XF(4) 

DS=DLQG((l.OD0+16,000*HHQ)/(i.OD0+8,0D0*RHO)) 

WRITEC6,13) 

NR1TE(o,15)  DS 
WRITE (6, 13) 

C1=DEXP(-(DIVL+DS) ) 

C2=2.0D0*C1 

RHS=DSdRT(4.0D04Cl*(1.0D0+Cl)) 

CALL  ASINH(X10,RHS) 

C4=QLOG(2.0DU)-OS 

C5sxl0*XI0 

XIB0=2.0D0*XI0 

C 

C  CALCULATE  XIBM  -  COORDINATE  OF  DOWNSTREAM  BOUNDARY  IN  XI  BAR  - 

C  ETA  BAR  PLANE, 

C 

XBE3XE+C4 

TERMsDEXP(XBt)-1.0D0 

RHSsUSORT(TERM*TERM-1.0D0) 

CALL  ASINH(X1E,RHS) 

XIBR=XIE4(1.0DO+C5/CPISOWXIEAXIE)) 

C 

CALL  SHEAR(RAD) 

KK  =  K 

CALL  XTGRIDCKK ,RAD) 

50  CONTINUE 
STOP 
END 

SUBROUTINE  SHEAR(RAD) 

C**«* ****** ********************************************************** 

C  THIS  SUBROUTINE  CALCULATES  SB AK  VS,  XI  BAR,  TO  BE  USED  IN  THE 

C  SHEARING  TRANSFORMATION. 

C*** ***************************************************************** 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  /BLKO 1 /  IMAX , JMAX , ITE , ITEM , ILAST , ISEG1 , ISEG2 
COMMON  /BLK02/  X I bM , XI 0 , XI BO 
COMMON  /BLK03/  Cl , C2 ,C3 ,C4 , C5 , PI , PISQ 
COMMON  /BLK04/  XU 101 ) , YF( 101 ) 

COMMON  /BLK05/  XIBC 1 51) , SBAR ( 1 51 ) , ETABLC 1 5 1 } 
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1230 

C 

1240 

DIMENSION  ETABU (151) 

1250 

C 

1260 

10 

FQRMATC5X, 'SHEARING  BOUNDARY  IN  XIBAR  -  ETABAR  PLANE*) 

1270 

11 

FORMAT ( 1  HO ) 

1280 

12 

F08MATC5X, 'I ' ,©X, 'XIBAR* »  9X , • ETA6L ' , 9X , ' ETA  BO ' , 9X , ' SBAR • ) 

1290 

13 

F0RMATCI6.4D14.4) 

1300 

14 

FORMAT(lriO,4X, 'UNABLE  TO  CONVERGE  XI  IN  50  ITER ATIONS ' /5X , 

1310 

1 

'  XIBAR  =»,D14.4) 

1320 

C 

1330 

C 

COMPUTE  NORMALIZED  AIRFOIL  COORDINATES  FOR  GIVEN  CYLINDRICAL 

1340 

C 

RADIUS  AND  TRANSFORM  TO  XI  BAR  -  ETA  BAR  PLANE.  THIS  STEP 

1350 

C 

GIVES  THE  FIRST  PONTIO.J  OF  THE  UPPER  BOUNDARY. 

1360 

C 

1370 

WRITE ( 6.11) 

1380 

WRITE ( 6.10) 

1390 

NRITEC6.il) 

1400 

T1=RAD*HAD 

1410 

DO  50  1=1, ITE 

1420 

YFl=YF ( 1 ) 

1430 

THF=DASIN(YFi/RAD) 

1440 

XBF I=XF ( I ) +  C4 

1450 

TriBFI=4.0D0*THF 

1 4b0 

C 

1470 

T2=0EXP ( X3F I ) 

1480 

P8AR=1.0DO-T2*DCCjS(THBFI) 

1490 

0BAR  =  T2*DSINCTHBH) 

1500 

OSU=UBAR*OHAR 

1510 

BETA=1.0D0-PBAR*PBAR-Q3(J 

1520 

ALPHA=DSGRT(BETA*EETA+4.0DO*QSO) 

1530 

RHSsDSORTCO. SCO* (ALPHA-BETA)) 

1540 

CALL  ASINH(XIF.RHS) 

1550 

ARG=PBAR/DC05H(XIF) 

1560 

IFC (ARG+l.ODO) .LT.O.ODO)  ARG=-1.0D0 

1570 

ETAF=DAC0S(ARG) 

1580 

C 

1590 

XMU=C5/(XIF*XIF+ETAF*ETAF) 

1600 

XIB(I)=XIF*(1.0DQ+XMU) 

1610 

ETABU(I)=ETAF»C1.0D0-XMU) 

1620 

50 

CONTINUE 

1630 

C 

1640 

C 

CONTINUE  UPPER  BOUNDARY  CALCULATION  BEYOND  AIRFOIL  T.E, 

1650 

c 

TO  XIBM . 

1660 

c 

1670 

DXIB=0.2D0 

1680 

ILAST=ITE+(XIBM-XIB(ITE) )/DXIB 

1690 

IFCILAST.GT.lSl)  ILAST=151 

1700 

ITEP=ITE+1 

1710 

WRITE ( 6 , 1 2 ) 

1720 

DO  100  lalTEP.ILAST 

1730 

XIBAR=XIb(I-l)+DXIB 

1740 

XIBCI )=XIBAR 

1750 

XXLsXIBAK 

1760 

DO  70  1T*1,50 

1770 

XMU*C5/(PIS0+XIL»XIL) 

1780 

XI*XlBAR/C1.0DO+XMU) 

1790 

IF(DAbS(XI-XIL).LT.1.0D-08)  GO  TO  80 

1800 

70 

XIL*XI 

1810 

WRITE (6, 14)  XIBAR 

1820 

STOP 

1830 

80 

ETABU(I)»PI*(1,0D0«XMU) 
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1840 

1850 

I860 

1870 

1880 
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1900 

1910 

1920 
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2060 
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2100 

2110 

2120 
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2160 
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2180 

2190 

2200 

2210 
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2250 

2260 

2270 

2280 
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2300 

2310 

2320 

2330 

2340 

2350 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 


100  CONTINUE 
C 

C  CALCULATE  LOWER  BOUNDARY  IN  XI  BAR  -  ETA  BAR  PLANE  AND 

C  SBAR. 

C 

DO  200  I=1,ILAST 
X IBAR=X IB ( 1 ) 

IF(XIBAR.GE.AIBO)  GO  TO  140 

XIL*XIBAR 

XMU  =  1 .000 

DO  120  11=1,50 

XI=XI6AP/(1.0D0*AMU) 

ARG=DCQSH(XI)-C2 

ETA=DACOS(AHG) 

XMU»C5/(XI*XI+ETA*ETA) 

IF(DA6S(XI-XIL).LT«1.0D-08)  GO  TO  130 
120  XIL=X1 

WRITE(6,14)  XIHAR 
STOP 

130  ETABL(I)=£TA*(1 .ODO-XMU) 

GO  TO  150 

140  ETA8L(I)=0.0D0 

150  5BAS(I)=ETABU(I)"LTABLCI) 

WRITE (b»13)  I,XI6(I),ETABL(I),ETABU(I),SBAR(I) 

200  CONTINUE 
RETURN 
END 

SUBROUTINE  X1GRID ( K , RAD ) 

C  THIS  SUBROUTINE  CALCULATES  THE  GRID  IN  THE  XBAK  -  THETABAR 

C  PLANE. 

C*9***9*******y*«999f***«****«********»*****¥***********«******** 

IMPLICIT  HEAL*8  (A-H,0-Z) 

COMMON  /BLKO 1/  I MAX , JK AX , HE , ITEM , ILAST , I SEG 1 , I SEG2 
COMMQi*  /BLK02/  XIBM , XIQ , XIBO 
COMMON  /BLK03/  C1,C2,C3,C4,C5,PI,PI5Q 
COMMON  /6LK05/  XI8C1S1),SBAR(151),ETABL(151) 

COMMON  /bLKOo/  SYO.SYl ,SZO,SXO,SSR 
COMMON  /BLK07 /  ZCC151) ,BIGZ(151) 

C 

DIMENSION  BIGX(lbl)  ,BIGU151  J 
DIMENSION  XCC151) ,YC1151) 

C 

11  F0RMATC1H0) 

12  FORMAT (5X, 'I' ,5X, 'J' ,SX, 'X' ,6X, 'XC' ,12X,'YC' , 1 2X , » ZC • , 1 2X , 
l'R',13X,,X,,13X,* THETA ' ) 

13  F0RMAT(3I6,6D14,4) 

C 

C  SET  UP  GRID  IN  COMPUTATIONAL  PLANE, 

C 

XBTF*XI3UIE) 

WPTE= ISEG1+I5EG2+ 1 

CALL  STRFX(XC,BIGX,ISEG1 ,1SEG2,IMAX,SX0#SSR,XIB0 , XBTE ,  X  IBM } 
CALL  5IRFY(YC,BIGY,JMAX,SY0,SY1) 

C 

C  DETERMINE  GRID  IN  PHYSICAL  PLANE. 

C 

1  =  1 

XIBAR=0,ODO 
S8I=SBAR ( 1 ) 

ETBLI*ETABL(1) 
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WRirE(6,U) 

wRlTE(o,12) 

WRITE ( o , 1 1 ) 

DO  70  J*1,JMAX 
ETA8AR=ETBLI+SBI*BIGY (J) 

P=U.25uO*£TAbAK*ETABAR+C5 

XIsO.OUO 

ETA=0.500*ETABAR+D5QKT(P) 

XBAR=DLOG(1.0DO-DCOS(ETA) ) 

THBAK=O.ODO 
XX  =  Xdftri-C4 
THETA=0 . 000 

WRITE (6, 13)  I,J,K,XC(I) , YC( J) ,ZCCK) , RAD , XX , THETA 
70  CONTINUE 
IBEG=1 
IEND=ITE 
DO  100  1*2  » IMAX 
XlBAwsblGXd) 

I F  C I . LE . N PTE )  GO  TO  80 
I8EG*ITE 
IEND=ILAST 
C 

C  INTERPOLATE  TO  FIND  5BAR  AND  ETABL  CORRESPONDING  TO  XIBAR. 

C 

80  CALL  If«TERP(XIB,SBAR,  XIBAR,  SBI , IBEG , IEND , INT , 0 ) 

CALL  iNTERPf XI B, ETABL, XIBAR, ETfaLI, IBEG, IEND, 1  NT, 1) 

C 

WRITEC6, 11) 

WRITE (6,12) 

WRITE ( 6,11) 

DO  100  J*1,JMAX 
ETABAR=ETBLI+SBI*BIGYCJ) 

Q*U.25U0*XIBAR*ETABAR 

P*0.25D0*(XIbAK*XIBAR-ETABAR*ETABAR)-C5 

XMU*D5QRT(P*P*4.0DQ*()*Q) 

XI=0.SD0*XiaAR+D6ORT(0.500*(XMU+P) ) 
ETA*0.5D0*ETABAR+DSORTC0.5D0*(XMU-P) ) 

C 

TlsOCOSH(Xl) 

T2=DCUS ( ETA ) 

ARG1*T1-T2 

XBAR=DLQG(ARG1) 

THBAR*DACOS((1,ODO-T1*T2)/ARG1) 

THETA*0.25D0*THBAR 

XX*X6AR-C4 

WRITE (6, 13)  I,J,K,XC(I) , YC ( J ) , ZC ( K ) , RAD, XX , THETA 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  ASINH ( ARG , RHS ) 

C***************** ************************************** ********* 

C  THIS  SUBROUTINE  COMPUTES  THE  INVERSE  HYPERBOLIC  SINE  USING 

C  NEWTON'S  METHOD. 

C»6**** *********************************************** *********** 
IMPLICIT  HEAL*8  (A-H.O-Z) 

C 

10  FOHMAT(1HO,4X, 'INVERSE  HYPERBOLIC  SINE  CALCULATION  FAILED  FOR 
lSINH(X)  *' , D1 4 . 7 ) 

C 

TE5T*DABS(RH5) 

IF (TEST .GT.l.ODO)  GO  TO  30 
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ARG=RHS 
GO  TO  40 

30  ARG=DLCGC2.QD0*TtST)*DSIGN(1.0D0,RHS) 

40  CONTINUE 

DO  50  *=1,50 
f A=DSInH(AKG)-RHS 
FPA=UCOSH(ARG) 

DAKG  =  -fr  A/FPA 

IF(DAbS(DARG) ,LT. 1 .00-10}  RETURN 
ARG=ARG+OARG 
SO  CONTINUE 

WRITE (6, 10)  KHS 

RETURN 

END 

SUBROUTINE  FOIL 

£***«*************♦************************************************ 

C  THIS  SUBROUTINE  GENERATES  (X,Y)  COORDINATES  FOR  A  SYMMETRIC 

C  JOUKOWSKY  AIRFOIL. 

q********** *************  ****************  ************  *************** 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON  /BLK01/  IMAX , JMAX , ITE, ITEM , ILAST, ISEG1 , ISEG2 
COMMON  /BLK03/  Cl , C2 , CJ , C4 , C5 , RI , PISU 
COMMON  /8LK04/  XF(lOi) , Yf (101) 

C 

10  FORMATCbX, 'AIRFOIL  COORDINATES') 

11  FORMAT  Cl  HO) 

12  FORMAT C 5X, ' I ' , 6X , • XF ' , 1 JX , ' YF ' ) 

13  FORMAT (16,2014.4) 

C 

DTH=PI/ITEM 
XFC1)=0.0D0 
YF(1)=0.0D0 
DO  50  1=2, ITEM 
Trfs ( 1-1 ) ♦DTH 
TlsDCOS(TH) 

XF(I)=0.5DQ*{1.0D0-T1) 

YF(I)=C3*(1.0DO+T1)*DSIN(TH) 

50  CONTINUE 

XF ( ITE ) =1 , 0D0 
YF ( ITE) =0 . ODO 
WRITE  (b, H) 

WRITE ( 6 , 1 0 ) 
w  RITE ( 0 , 1 1 ) 

»RIT£ (6,12) 

DO  60  1=1, ITE 

60  WRITE (6, 13)  1,XF(I) ,YF(I) 

RETURN 

END 

SUBROUTINE  INTERP ( XX , Y Y , X INT , Y I  NT , I BEG , I END , I  NT , ISW ) 

C********* ***************  4  *********************  **************** 

C  THIS  SUBROUTINE  USES  LAGRANGE  CUBIC  INTERPOLATION  TO 

C  DETERMINE  YINT  FOR  A  GIVEN  XINT. 

C 

c  xx  =  independent  variable. 

C  YY  =  DEPENDENT  VARIABLE. 

C  1BEG  =  INITIAL  INDEX  FOR  INTERPOLATION  RANGE, 

C  I END  =  FINAL  INDEX  FOR  INTERPOLATION  RANGE, 

C  INT  =  UPPER  INDEX  OF  INTERPOLATION  INTERVAL, 

C  IS*.  =  INTERPOLATION  INTERVAL  SEARCH  SWITCH, 

C  0  PERFORM  SEARCH, 

C  1  OMIT  SEARCH. 
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c******«**********»*****************************«***«**********« 

IMPLICIT  HEAL*8  (A-H,G-Z) 

C 

DIMENSION  XX (151), YY (151) 

C 

IF(ISw.GT.O)  GU  TO  75 
60  00  70  I=1SEG,IEND 
INT  =  I 

IF(XXU).GT.XINT)  GO  TO  75 
70  CONTINUE 

75  IF(lNT.EQ.(IbEG+l) )  GO  TO  80 
IF(INT.EG.IEND)  GO  TO  90 
Il=INT-2 
1 2  =  1  N  T— 1 
I3=INT 
I4=INT+1 
GO  TO  100 
90  I1=IBEG 
I2=IbEG+ 1 
I3=IBEG+2 
I4&IBEG+3 
GO  TO  100 
90  I1=I£ND«3 
12= I EnD«2 
I3=IEnO-1 
1 4  =  I  £  N  D 
loo  continue 
X1=XX(I1) 

X2=XX ( 12) 

X3=XX ( 13) 

X4=XX(I4) 

CFl=(XIwT-X2)*(XlNT-X3)*(XINT-X4)/((Xl-X2)*(Xl-X3)*(Xl-X4) ) 
CF2=(XINT-X1)*(XINT-X3)*(XINT-X4)/((X2-X1)4(X2-X3)*(X2-X4) ) 
CF3*(XI.Vf-Xl)*(XlNT-X2)*(XlMT-X4)/((X3-Xl)*(X3-X2)*(X3«X4) ) 
CF4=(XINT-Xl)*(XlNT-X2)»(XIM-X3)/(  ( X9-X1 ) * ( X4-X2 ) * ( X4-X3 )  ) 
YINT=CFl4YY(Il)+CF21»YY(I2)+CFi4YY(I3)+CF4*YY(I4) 

RETURN 

END 

SUBROUTINE  STHFY ( X I , T , NPT , S YO , 5 Y 1 ) 

C  THIS  SUBROUTINE  GENERATES  A  NONUNIFORM  POINT  DISTRIBUTION 

C  USING  VINOKURS  TWO-SIDED  STRETCHING  FUNCTION,  AS  GIVEN  IN 

C  NASA  CR-3133, 

C*******¥***«*********************************************f***f** 

IMPLICIT  REAL*8  (A-H,0-Z) 

C 

DIMENSION  XI ( 1 51 ) # T ( 1 5 1 ) 

C 

C  COMPUTE  XI. 

C 

DXI*1.0DQ/(NPT-1) 

DO  40  J=1 , NPT 
40  XI(U)=(J-1)*0XI 
C 

C  COMPUTE  DELTA  Y. 

C 

AsDSQRT ( SY0/SY1 ) 

B»DSQRT(SY0*SY1 ) 

TEST* 2. 7 8296 8  IDO 
IF(B.GT.TEST)  GO  TO  50 
Y8AR=B-1.0D0 


-51- 


4280 

4290 

4300 

4310 

4320 

4330 

4340 

4350 

4360 

4370 

4300 

4390 

4400 

4410 

4420 

4430 

4440 

4450 

4460 

4470 

4400 

4490 

4500 

4510 

4520 

4530 

4540 

4550 

4560 

4570 

4500 

4590 

4600 

4610 

4620 

4630 

4640 

4650 

4660 

4670 

4680 

4690 

4700 

4710 

4720 

4730 

4740 

4750 

4760 

4770 

4780 

4790 

4000 

4010 

4820 

4030 

4840 

4850 

4860 

4070 

4080 


30  March  1983 
GHH : lhm 

DELY=(  (  (  (C-0. 00 107941 23U0*YBAR+0. 007 7424461  DO)  *Y BAR 
1-0.024907295DO)*YBAR+0.057321429DO)*YBAR-0. 15D0)*YBAR 
2+1.0D0)*DSURT(6.0D0*YSAR) 

GO  TO  6u 
50  V=DL0G(b) 

0*1 .0D0/8-0. 02  8 527  4  3  IDO 

DELY=( ( (8. 5679591 1 D04*-2. 6294547  DO )*0  +  l .949644 JD0)*H 
1 +0.24902722D0)*w-0. 02041 7 93D0+V+C1 .000+1 . QDO/VJ* 

2DL0G ( 2 .0004V ) 

60  CONTINUE 
C 

C  COMPUTE  T. 

C 

C1=A*DSINH(DELY) 

C2=1.0D0-A*0COSHCDELY) 

DO  70  1*1, NPT 
''  =  DTANH(DELY*XI(I)) 

T(I)eFN/(Cl+t2*FN) 

70  CONTINUE 
RETURN 
END 

SUBROUTINE  STRKX(XI,T,NSEG1,NSEG2,NMAX,SX0,SSR,XIB0,XBTE,XIBM) 
C*********f*************9***t**>M*****************t***V*«*********** 

C  THIS  SUBROUTINE  GENERATES  A  NONUNIFOHM  POINT  DISTRIBUTION 

C  SPECIALIZED  TO  THE  COORDINATE  TRAPPED  AROUnU  THE  AIRFOIL. 

£***********4 i********9*********************v****ff****¥**«?********* 

IMPLICIT  REAL*b  (A-H,0-Z) 

C 

DIMENSION  X I ( 1 5 1 ) , T  C 1 5  1 ) 

C 

C  SEGMENT  NUMBER  1. 

C 

TTEsXBTE/XIBO 

TMAX=XIBM/XIB0 

dxi=i.odo/nsegi 

NP1=NSEG1+1 
S1=0.5D0*(SX0-1 .ODO) 

DO  50  1=1, NP1 
XX=U-1  )*DXI 
XI(I)=XX 

50  T(I)=XX*(1.0D0  +  S1M1 . 0D0-XX3  * ( 2  .  ODO-XX )  ) 

c 

C  SEGMENT  NUMBER  2. 

C 

AA*0. 5004(3. ODO-SXO) 

XWTE=NSEG2*DAI 

BB*(TTE-1 ,0DO-AA4X«TE)/(XWTE4X>tTE) 

NP2*NSEG2+1 
DO  60  K=2 , NP2 
I *NSEG 1 +K 
XW*(A-1 )4DXI 
XI ( I ) S1 , ODO+X0 
60  T(I)*1.0DO4XW*(AA  +  X'4*BB) 
c 

C  SEGMENT  NUMBER  3. 

C 

N3*N5EG1+N5EG2 
NP3*N3+1 
X ITE*XI ( NP3 ) 

DT1*T(NP3)-T(N3) 

S1*OT1/(SSH-1,ODO) 


GHH :  lhm 


4690 

KMAX*151-NP3 

4900 

DO  70  K=2,KMAX 

4910 

ISN3+K 

4920 

XI(I)sXITE+(K-l)*DXI 

4930 

ti»tte+si*cssr**(k-i)-i .ono) 

4940 

T(I)=TI 

4950 

IF(TI.GE.THAX)  GO  10  80 

4960 

70 

CONTINUE 

4970 

80 

NMAX  =  i 

4960 

C 

4990 

C 

RESCALE  VARIAbLES. 

5000 

C 

5010 

SCALE=XtlTE/XITE 

5020 

DO  90  I=1,NMAX 

5030 

XI(I)=SCALE*XI(I) 

5040 

90 

T(I)=XIflO*T(i) 

5050 

RETURN 

5060 

END 

5070 

SUBROUTINE  S1’RFZ(XI,T,NPT,S0) 

5080 

C**** 

***#******»************«*¥***»***V****************************** 

5090 

c 

THIS  SUBROUTINE  GENERATES  A  NQNUNIKOKM  POINT  DISTRIBUTION  USING 

5100 

c 

VINOKURS  ONE-SIDED  STRETCHING  FUNCTION. 

5110 

c**»* 

5120 

IMPLICIT  REAL*8  (A-H,0-Z) 

5130 

c 

5140 

DIMENSION  XK151)  ,T(151) 

5150 

c 

5160 

c 

COMPUTE  XI. 

5170 

c 

5180 

DXIsl.ODO/CNPT-l) 

5190 

DO  40  K=1,NPT 

5200 

40 

XI(K)a(K-l)*DXI 

5210 

c 

5220 

c 

COMPUTE  DELTA  Y. 

5230 

c 

5240 

TEST=2.7829661D0 

5250 

IF(SO.GT.TEST)  GO  TO  50 

5260 

YBAR=50-1.0DO 

5270 

DELY=( ( (  C (-0.0010794 l23U0*¥BAR  +  0. 007 7424461 DO)* YBAR 

5280 

1-0. 024907 295D0)*IPAR+0. 057 321 429 DO ) * YBAR-0 , 1 5D0 ) * YBAR 

5290 

2+1. 0D0) *DSUR  T(6.0D0*YBAR) 

5300 

GO  TO  60 

5310 

50 

VsDLOG(SO) 

5320 

W  =  1 .0D0/S0-0. 028527 43  IDO 

5330 

DELY*( ( (8. 5679591 1 DO+W-2 . 6294547D0 ) *W+ 1 . 9496443D0 ) *M 

5340 

1+0,24902722DO)*R-0.02U41793DO+V+(1 .ODG+1 .ODO/V)* 

5350 

2DLOG12,ODO*V) 

5360 

60 

CONTINUE 

5370 

c 

5300 

c 

COMPUTE  r. 

5390 

c 

5400 

C1=0.5D0*DELY 

5410 

C2=l.QD0/D'fAUH(Cl) 

5420 

DO  70  Ksl.NPT 

5430 

T(K)sl.0D0+C2*DTANH(Cl*(XI(K)-1.0D0)) 

5440 

70 

CONTINUE 

5450 

RETURN 

5460 

END 
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